# STAT 415/615  Regression (M. Baron)

# Python Lab 4. Inference for the Univariate Linear Regression

# Let’s work with textbook data on Copier Maintenance, problems 1.20, 2.5, etc. The data are available on our Blackboard, but also on the internet in the public domain. Here is how we can read them in one step.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats

# Read the data directly from the internet:

url = "http://www.cnachtsheim-text.csom.umn.edu/Kutner/Chapter%20%201%20Data%20Sets/CH01PR20.txt"

C = pd.read_csv(
    url,
    sep=r"\s+",
    header=None,
    names=["Y", "X"]
)

# Look at the first few observations:

C.head()

# Assign X and Y

X = C["X"]
Y = C["Y"]

# Fit a regression model of Y on X:

reg = smf.ols("Y ~ X", data=C).fit()

# Summary of regression inference

# The regression summary includes t-tests for the regression coefficients, analysis of residuals, the F-test, and R².

print(reg.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      Y   R-squared:                       0.957
Model:                            OLS   Adj. R-squared:                  0.957
Method:                 Least Squares   F-statistic:                     968.7
Date:                Sun, 30 Aug 2026   Prob (F-statistic):           4.01e-31
Time:                        00:03:13   Log-Likelihood:                -161.27
No. Observations:                  45   AIC:                             326.5
Df Residuals:                      43   BIC:                             330.2
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -0.5802      2.804     -0.207      0.837      -6.235       5.075
X             15.0352      0.483     31.123      0.000      14.061      16.009
==============================================================================
Omnibus:                        1.399   Durbin-Watson:                   2.402
Prob(Omnibus):                  0.497   Jarque-Bera (JB):                1.167
Skew:                          -0.388   Prob(JB):                        0.558
Kurtosis:                       2.853   Cond. No.                         12.5
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

# Inference about regression coefficients

# The fitted regression model contains several components. We can inspect the estimated coefficients directly:

reg.params

Intercept    -0.580157
X            15.035248
dtype: float64

# We can save the intercept and slope:

intercept = reg.params["Intercept"]
slope = reg.params["X"]

slope

15.035248041775459

# Confidence intervals for regression coefficients

# By default, `conf_int()` gives 95% confidence intervals for β₀ and β₁:

reg.conf_int()

# We can set the desired confidence level.

# A 90% confidence interval:

reg.conf_int(alpha=0.10)

# A 99% confidence interval would be wider:

reg.conf_int(alpha=0.01)

# Estimation of mean responses and prediction of individual responses

# Predict the response Y for X = 3:

new_data = pd.DataFrame({"X": [3]})

reg.predict(new_data)

0    44.525587
dtype: float64

# Confidence interval for E(Y) at X = 3:

prediction = reg.get_prediction(new_data)
prediction.summary_frame(alpha=0.05)

# The resulting table contains both the predicted mean response and the confidence interval for the mean response.

# For a prediction interval for an individual future observation, use the same prediction results:

prediction.summary_frame(alpha=0.05)[
    ["mean", "mean_ci_lower", "mean_ci_upper",
     "obs_ci_lower", "obs_ci_upper"]
]

# Here:

#  `mean` is the predicted value
#  `mean_ci_lower` and `mean_ci_upper` give the confidence interval for E(Y) at X = 3
#  `obs_ci_lower` and `obs_ci_upper` give the prediction interval for an individual Y at X = 3

# Confidence band for the whole regression line

# We can construct a confidence band for the entire regression line using the Working–Hotelling method.

# First, obtain the sample size:

n = len(Y)

# Save the residuals:

e = reg.resid

# The estimated standard deviation is the square root of the mean squared error:

s = np.sqrt(np.sum(e**2) / (n - 2))

# Calculate the Working–Hotelling multiplier:

W = np.sqrt(2 * stats.f.ppf(0.95, 2, n - 2))

# Obtain the fitted values:

Yhat = reg.fittedvalues

# Calculate Sxx:

Sxx = (n - 1) * X.var()

# Calculate the upper and lower confidence bands:

upper_band = (
    Yhat
    + W * s * np.sqrt(
        1/n + (X - X.mean())**2 / Sxx
    )
)

lower_band = (
    Yhat
    - W * s * np.sqrt(
        1/n + (X - X.mean())**2 / Sxx
    )
)

# Plot the observations and the regression line:

plt.scatter(X, Y)
plt.plot(X, Yhat, linewidth=2)

plt.xlabel("X = number of copiers")
plt.ylabel("Y = service time")

plt.show()

# Add the confidence bands:

plt.scatter(X, Y)
plt.plot(X, Yhat, linewidth=2)
plt.plot(X, upper_band, linewidth=1)
plt.plot(X, lower_band, linewidth=1)

plt.xlabel("X = number of copiers")
plt.ylabel("Y = service time")

plt.show()

# Regression residuals

# The regression residuals are the differences between the observed and fitted values:

e = reg.resid

# We can save the residuals for further analysis.

# A summary of the residuals can be obtained with:

e.describe()

# Plot the residuals:

plt.plot(e, "o")
plt.axhline(0, linewidth=1)
plt.xlabel("Observation")
plt.ylabel("Residual")

plt.show()

# We can verify that the sum of the residuals is zero (up to numerical rounding):

e.sum()

-1.1013412404281553e-13

# The residuals are also orthogonal to X. Thus, their cross-product with X is zero (up to numerical rounding):

np.sum(e * X)

-1.1368683772161603e-12

# ANOVA: sums of squares and F-test

# The ANOVA table for the regression can be obtained using:

sm.stats.anova_lm(reg)

# The ANOVA table contains:

# Degrees of freedom (`df`)
# Sum of squares (`sum_sq`)
# Mean squares (`mean_sq`)
# F-statistic (`F`)
# p-value (`PR(>F)`)

# For a univariate linear regression, the ANOVA F-test tests the null hypothesis H₀: β₁ = 0
# against the alternative Hₐ: β₁ ≠ 0.

# The F-test and the t-test for the slope are equivalent when there is only one predictor. In fact, F = t².

# Thus, the ANOVA provides another way to test whether the predictor X is significantly associated with the response Y.

# What does NaN mean here? For the Residual row, the F-statistic and p-value are simply not applicable.

	0	1
Intercept	-6.234843	5.074529
X	14.061010	16.009486

	0	1
Intercept	-5.293780	4.133467
X	14.223144	15.847352

	0	1
Intercept	-8.137064	6.976751
X	13.733279	16.337217

	df	sum_sq	mean_sq	F	PR(>F)
X	1.0	76960.422977	76960.422977	968.657196	4.009032e-31
Residual	43.0	3416.377023	79.450628	NaN	NaN